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PREPACE 

Part  of  the  Project  RAND  research  program  consists  of 
basic  supporting  studies  In  mathematics.  One-  aspect  of  this 
Is  concerned  with  optimization  processes.  In  this  field,  the 
technique  of  dynamic  programming  has  developed  Into  a  powerful 
mathematical  tool. 

In  the  present  Memorandum  the  authors  Investigate  the 
applicability  of  polynomial  approximation  as  an  adjunct  to 
the  technique  of  dynamic  programming. 
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SUMMARY 

In  this  Memorandum  the  authors  initiate  the  application 
of  the  simple  yet  powerful  computational  technique  of 
polynomial  approAiiriation  to  problems  in  dynamic  programming. 

The  theoretical  applicability  of  orthogonal  polynomials 
is  first  discussed  and  then  applied  to  one—  and  two-dimensional 
allocation  problems.  Numerical  results  obtained  from  FORTPIAN 
programs  involving  Legendre  polynomials  are  presented. 
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POLYNOMIAL  APPROXIMATION — A  NEW  COMPUTATIONAL  TECHNIQUE 
IN  DYNAMIC  PROGRAMMING — I:  ALLOCATION  PROCESSES 


I ■  INTRODUCTION 

The  problem  oi’  maximizing  the  function 

(1.1)  F(x^,a2,  .  .  +  SaIx;,)  +  •••  + 

over  the  domain 


(1.2)  "t  +  *  *  *  +  i  ^ 


can  be  reduced  via  familiar  dynamic  programming  techniques 
(see  [ij)  to  that  of  determining  the  sequence  of  functions 
(f^(A)]  generated  by  means  of  the  recurrence  relation 


(l.Q)  f, (x)  =  g, (x), 


fn^^(x)  =  max  [g„^T(y)  +  fjx  -  y)  j 


0<y<x 


n+1' 


The  problem  can  thus  oe  solved  numerically  in  a  very  simple 
fashion,  regardless  of  the  complexity  of  the  functions  g^(x). 
A  number  of  Lrr.portant  allocation  processes  can  be  resolved  in 
this  way.  If  we  consider  cases  in  which  two  distinct  types  of 
resources  must  be  allocated,  we  face  the  problem  of  maximizing 
the  function 

F(x^,x.,,  .  .  .,Xj^;y^,y2,  .  .  .,yj^) 

=  +  g2(x2..yJ  +  •••  +  SN(^^N^yN^ 


(1.4) 
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over  the  domain 

(1.5)  +  •  •  ♦  +  =  X,  Xj^  ^  0, 

Yl  +  y2  +  • • •  +  Yn  =  Yi  2  0. 

TheoretlcallY,  there  is  no  difficultY  in  applying  the  same 
techniques . 

The  maximization  problem  can  be  reduced  to  the  deter^ 
mination  of  the  sequence  {f^(x,y)]  generated  by  means  of 
the  recurrence  relation 

(1-^')  ^  ^n^^  -  w,y  -  r)]. 

CK  w<x 

In  pi’inciple,  this  equation  can  be  solved  computationally 
using  the  same  technique  that  applies  so  well  to  (I.5),  In 
practice  (see  [l]  for  a  discussion),  questions  of  tiJTie  and 
accuracy  arise.  There  are  a  number  of  ways  of  circumventing 
these  difficulties,  among  which  the  Lagrange  multiplier  plays 
a  significant  role. 

In  this  series  of  papers,  we  wish  to  present  a  number  of 
applications  of  a  new,  simple  and  quite  powerful  method,  that 
of  polynomial  approximation.  We  shall  begin  with  a  discussion 
ol'  the  allocation  process  posed  in  the  foregoing  paragraphs 
and  continue,  in  subsequent  papers,  with  a  treatment  of 
realistic  trajectory  and  guidance  processes.  In  a  separate 
series  of  papers  we  shall  apply  this  fundamental  attack  upon 
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dimensionallty  to  the  solution  of  a  number  of  the  equations 
of  mathematical  physics. 

We  would  like  to  express  our  appreciation  to  Oliver  Gross 
for  the  analytic  solution  of  some  test  problems  we  used  to 
check  the  accuracy  of  our  techniques,  and  for  his  general 
interest  in  the  program. 

_ POLYNOMIAL  APPROXIMATION 

In  the  systematic  study  of  dynaiTiic  programming  as  a 
computational  algorithm  given  in  [l],a  function  f(x)  is 
considered  to  be  a  table  of  values  at  an  appropriate  set  of 
grid  points: 


X 

0 

0 

1 

A 

1 

2A 

^2 

KA  1 

In  this  way,  the  function  f{x)  is  stored  in  the  computer. 

For  functions  of  one  variable  with  K  =  100  or  1000,  this 
is  a  reasonably  convenient  way  to  proceed.  For  functions  of 
two  variables,  this  procedure  becomes  a  bit  Inconvenient  since 
(K  +  1)  values  for  x  combined  with  (K  +  l)  values  for  y 
yields  a  total  set  of  (K  +  1)“^  values.  Consequently,  when 


we  encounter  functions  of  three  or  more  variables,  we  must 
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balance  accuracy  against  time  and  the  limited  storage  of 
contemporary  computers  in  our  choice  of  K. 

The  storage  of  functional  values  by  means  of  a  table  of 
values  is  ideally  suited  to  the  treatment  of  problems  involving 
functions  of  quite  arbitrary  form.  It  is,  on  the  other  hand, 
quite  wasteful  and  inefficient  if  we  are  dealing  with  functions 
possessing  a  definite  structure,  which  is  to  say,  situations 
in  which  there  is  a  high  correlation  between  the  values  of 
f(rA)  and  f(sA).  Since  functions  of  this  type  occur  in 
many  important  applications,  and  throughout  mathematical 
physics,  it  is  worthwhile  to  develop  methods  which  take 
advantage  of  the  "smoothness"  of  the  function. 

One  such  method  is  polynomial  approximation,  or  to  be 
precise,  generalized  polynomial  approximation.  We  represent 
the  function  in  the  form 

M 

(r.S)  f(x)  =  Z  cpi^(x), 

1 

where  the  are  known  elementary  functions  such  as  x^, 

sin  kx,  Pl,(x)  (the  Legendre  polynomial  of  degree  k)  or 

(the  Cebycev  polynomial  of  degree  k),  and  then  store 
the  function  for  all  values  of  Interest  by  means  of  the  set 
of  K  coefficients  [a^, a^,  • .  • , ] . 

It  is  important  to  point  out  that  (2.2)  is  particularly 
useful  in  automatically  furnishing  the  interpolation-values 
frequently  needed  in  dynamic  pT’ogramming  calculations.  If 


one  uses  a  table  of  values  of  the  forai  shown  in  (2.1), 
interpolation  is  frequently  a  source  of  difficulty. 

To  determine  the  coefficients  a^^  it  is  convenient  to 
make  the  be  an  orthonormal  set.  Then 

(2.3)  f  f(x)'Pj^(x)dx, 

0 

assuming  for  the  purposes  of  convenience  that  0  ^  x  <  1.  We 
can,  of  course,  use  a  Cebycev  fit  instead,  and  we  will  explore 
this  in  subsequent  papers.  A  priori,  we  would  suspect  that  it 
would  be  more  efficient  to  use  a  mean— square  fit  (implied  by 

(2.3) )j  and  take  M  to  be  larger,  than  to  go  to  the  trouble 

of  determining  the  to  minimize  the  function 

M 

(2.4)  max  !f(x)  -  Z  a,9-(x)i, 

0<x<l  k=l  ^ 

which  for  a  fixed  M  may  be  expected  to  yield  a  more  accurate 
approximation.  Alternatively,  one  could  use  an  approximation 
by  polygonal  functions  [l]. 

To  evaluate  the  without  requiring  a  knowledge  ol'  too 

many  values  of  f(x),  we  use  a  quadrature  technique 

R 

(2„5)  =  Z  w^f(x^)cp^(x^), 

i=il 

where  the  weights  w^,  and  the  quadrature  points  x^^  are 


chosen  so  that  the  equation 


(2.6) 


1  R  _ 

g(x)dx  =  Z  w.g(x, ) 
i=l  ^ 

is  exact  for  polynomials  in  x  of  degree  (2R  —  1)  or  less. 
This  requirement  narrows  us  down  to  the  Gauss  quadrature 
technique  [2],  If  we  use  generalized  polynomials 
(expressions  such  as  (2.b)),  we  will  obtain  different  weights 
and  quadrature  points. 

We  may  then  store  the  function  f(x)  for  0  <  x  <  1 
by  storing  the  coefficients  a^^  or  the  particular  values 
f(x^)  which  enable  us  to  compute  the  a,^. 

3.  APPLICATION  TO  DYN.4MIC  PR0GRAJV1.MING 

Consider  now  the  application  of  these  ideas  to  the 
computational  solution  of  the  functional  equations  of  dynamic 
programming.  Suppose  that  we  wish  to  compute  the  sequence 
(f^(x)}  determined  by 

(^•i)  [Sj^(y)  +  f^l(x  -  y)  ]., 

0^y<x 

N  2  given  that  f2(.!<-)  =  g]^(^)-  To  avoid  the  tabulation  of 
each  of  the  functions  x— grid  [O,  A,  .  .  . ,  KA  ]  j 

where  K  may  be  a  large  number,  we  approximate  to  each 
function  Tj^(x)  in  the  manner  indicated  above.  Starting 
with  f-j^(x),  we  store  the  values  ^  ~  1^2,..., R, 

needed  to  evaluate  T^(x  —  y)  in  the  formiula  determining 

(3.2)  f-(x)  -  max  [g.(x  )  +  f^(x-y)]. 

0<y<x 
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We  then  determine  successively 

a  set  of  values  which  stores  the  function  f2(x).  This 
process  is  then  repeated. 

Although  the  calculation  of  fj^_^(x  -  y)  using  (2.2) 
is  far  more  time-consuming  than  taking  a  value  from  storage, 
we  expect  to  gain  time  over— all  because  of  the  fact  that  we 
are  required  to  calculate  only  a  few  values,  f^(x^), 

1  =  l,2,...jR,  at  each  stage. 

i-_  THE  LEGENDRE  POLYNOMIALS 

Since  in  allocation  processes  we  have  a  range  [0,Xq] 
which  stays  constant  as  the  process  continues  from  stage  to 
stage,  we  can  normalize  a- d  consider  the  basic  Interval  to 
be  [0,l].  Since  the  interval  is  finite  and  we  want,  at  the 
moment,  a  polynomial  approximation,  we  shall  employ  Legendre 
polynomials . 

Let  f  (x)  denote  the  standard  Legendre  polynomial, 
defined  over  [—  1,1],  and  let  <Pj^(x)  be  defined  by  the 
relation 

(4.1)  ^j^(x)  =  (2k  +  l)^/2pj^(2x  -  1). 

The  seque^nce  (<Pj^(x))  is  then  orthonormal  over  [0,l],  i.e., 

r  1 

(4.2)  j  (Pj^(x)cp.(x)dx  = 

0 

From  the  standard  recurrence  relations  for  obtain 


the  relation 
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(4.3)  =  3^^^(2x  -  1), 


-  (21  + 

-  (i  +1 


(21  +  .1)^/2^2x  -  l)q)i+i(x) 


- TTTp-  ‘Pi(^') 

(21  -  1)-' "  ^ 


This  relation  makes  the  evaluation  of  the  sequence  of  values 
of  9j_(x)  for  any  x  a  relatively  simple  matter. 

EXWL^ 

In  order  to  experiment  with  this  new  approximating; 
procedure,  we  devised  a  FORTRAN  program  for  the  IFM— 7090 
whereby  at  each  stage,  after  obtaining  the  new  coefficients 
we  computed  f^(x)  from  (2.2)  for  a  succession  of 
as  many  values  of  x  as  desired,  and  printed  the  result  after 
each  computation.  We  used  two  modes  of  output,  either  a  list 
of  numerical  values,  [x,f^(x)],  or  a  graphical  plot  (done 
directly  by  the  7090)  of  x  versus  f^(x). 

We  experimented  with  several  types  of  functions  gj^(x) 

for  which  the  results  could  be  derived  from  analytic 

♦ 

considerations.  Using  the  known  analytic  results  as  a 
checkpoint,  we  varied  the  parameters  R,  M,  and  the  grid  size 
H,  in  order  to  determine  the  degree  of  accuracy  we  might 
expect  in  general.  We  found  remiarkably  good  agreement  to  2 


* 

Oliver  Gross  was  of  considerable  assistance  to  us  in 
this  respect. 
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or  more  significant  figures  with  a  relatively  low  order  of 
approximation,  namely  r  =  5,  m  =  6.  This  is  encouraging 
from  the  point  of  view  of  extending  the  method  in  the  experi¬ 
mental  investigations  of  higher-dimensional  allocation 
processes  where,  as  pointed  out  above,  time  and  storage 
aspects  become  significant.  We  also  Incorporated  in  our 
program  restraints  of  the  form  0  <  <  x,  since 

constraints  of  this  nature  often  occur  in  applications. 

a.  Time  Est.imates 

Following  are  some  estimates  of  the  execution  time 
required  on  the  7090: 

1.  10  seconds  for  4  stages,  R  =  9,  M  =  6  and  a  grid 

of  O.Od  both  for  the  search  and  output  listing. 

.  3  minutes  for  a  total  of  4  cases  of  10  stages  each; 

R  =  3,  M  =  4;  R  =  0,  M  =  6;  R  =  7,  M  =  8;  R  =  10,  M  =  11. 

The  grid  in  the  search  for  the  maximum  is  0.01  and  the  output 
listing  is  given  for  every  x  along  the  grid. 

3.  3  minutes  for  the  total  of  4  cases  mentioned  above, 
10  stages  per  case  and  a  grid  of  0.01  for  the  search.  The 
output  is  a  graph  where  the  Independent  variable  x  is 
listed  at  intervals  of  0.029- 

Numerical  Results 

1  '7' 

1*  =  i(-^)  ' ''  '  Using  the  Schwarz  Inequality,  we 

readily  obtain  the  values 

=  (f(N  +  1)(2N  +  l)x)^/^. 


For  R  ■  10,  M  »  11,  H  -  O.Ol  we  found  poor  agreement  at 
the  origin.  This  Is  to  be  expected  because  of  the  Infinite 
slope  at  X  »  0.  Agreement  between  exact  and  computed  values 
was  good  as  soon  as  x  moved  away  from  the  singular  value  0, 
as  can  be  seen  from  the  following  table: 


Function 

Exact 

Computed 

fi(o) 

0.0 

0.064 

fl(2) 

0.448 

0.447 

fl(l) 

1.00 

1.00 

fio«)) 

0.0 

3.13 

fio'*) 

19.6 

19.6 

=  l(x  +  Otl)  We  avoided  the  previous 

difficulty  at  x  ■  0  (see  the  computed  value  of 
above),  but  found  the  function  still  rather  sensitive  near 
the  origin.  As  N  (the  number  of  stages)  Increased,  the 
agreement  at  x  ^  0  decreased. 

3*  gj^(*)  ■  +  1)^'^^.  The  Schwarz  Inequality  yields 

the  upper  bound 

i  +  1)<2N  +  l)(x  +  N))^/^. 

However,  since  the  x^  are  subject  to  the  rastralnt 
0  ^  x^  ^  X,  we  do  not  necessarily  achieve  the  upper  bound. 
Some  exact  values  based  on  an  analysis  by  0.  A.  Gross  are 
listed  as  check  points,  R  >  10,  M  >  11,  H  «  0.01. 
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Function 

UDoer  bound 

Exact  value 

Computed  value 

fi(o) 

1.00 

1.00 

1.00 

fl(l) 

1.41 

1.41 

1.41 

f3(0) 

6.48 

6.00 

fjC-^) 

7.00 

6.67 

6.67 

f7(.55) 

52.1 

29.1 

29.1 

^10^^^ 

— 

59.5 

59.3 

4.  gi(x)  -  R  -  10.  W  -  11.  H  -  .01. 


Function 

Exact 

Computed 

0.196 

0.197 

f2(*7) 

0.659 

0.659 

0.203 

0.204 

f3(.9) 

0 .860 

0.862 

0.218 

0.217 

^10^^^ 

1.001 

1.001 

6.  two-dimensional  approximation 

'Hie  problem  of  maximizing  the  function 


(6.1)  g^Cx^.y^) 

over  the  domain 
N 

(6.2)  ^2^  x^  -  X.  0  1  a^ 

N 

2  -  y.  0  ^ 

1-1 


1  *1  <  ^1' 
S  yi  1  di» 


i  X. 

yi  i  y* 
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can  be  readily  handled  by  the  methods  described  in  Secs.  2 
and  5  • 

The  dynamic  programming  recurrence  relation  is: 


=  max 
R 


yj  + 


Vi<^  -  ^N-y  -  yw) 


where  R  is  determined  by  a^^,  <  <  min(x,bj^), 

Gv,  <  Yv,  <  min(y,d.,).  Let  x,,  i  ^  1,...,R,  be  the  roots 
N  —  N  —  i'<  1 

of  the  normalized  Legendre  polynomial  j2^„(x),  and  let  (yj 

^  J 

be  a  duplicate  set  of  these  quadrature  points.  Each  function 
fj^(x,y)  is  expressed  approximately  by  the  relation 

MM 

(6.3)  f^,(x,y)  =2  S  a  ^  ^  d  (x)  ^ Ay), 

r-1  s=l  ' 


where  the  coefficients,  using  the  quadrature  method,  are 
given  by 


(6.4) 


(N 


R  R 

ifi  jfi 


As  in  the  one— dimensional  case,  we  start  with  the  known 

values  f^(x,y)  =  gj^(x^,y^),  x^  =  x,  y^  =  y.  In  stage  n 

we  store  the  values  ^  and 

then  compute  and  store  the  values  a  in  the  storage 

r ;  s 

allotted  to  the  previous  stage.  The  latter  coefficients  are 
utilized  in  the  computation  of  ^’j^(-^  “  -^n+l^^  ”  ^n+1^ 
obtain  the  values  of  ^n+l^'^i’^J^  next  stage. 


-13- 


_ EXAMPLES 

a.  Time  Estimate 

The  execution  time  required  on  the  7090  was  2  minutes  for 
•'I  stages j  with  R  =  9,  M  =  fcj  a  grid  of  H  =  0.09  in  the  two- 
dimensional  search,  and  an  output  listing  of  3  test  values  of 
f(x,y)  in  each  stage. 

b.  Numerical  Results 

1.  g^(x,y)  =  (21  -  l)^/'‘'-(xy)^^^  •  Using  the  Holder 

V  ‘  1  ^’4 

inequality,  we  readily  obtain  the  exact  values  £^{x.,y)  =  n(xy)  ^ 
x^(x,y)  =  (^n-l)x/n‘\  yj.^(x,y)  =  (2n-l)y/n^  .  Oux'  results  for 
R  =  9^  M  =  b,  and  a  grid  of  0.05  in  the  search  are: 


Function  j 

I  Exact 

Computed 

f2(o,.5) 

1.40 

1 .40 

[ 

0075 

0O5  (to  the 

9(1*1) 

4.00  1 

5.91 

nearest  0*09) 


Here  again  the  agreement  was  poor  at  the  origin,  presumably 
because  of  the  singular  behavior  of  (xy)^'^^  at  x  =  0,  y  =  0. 
To  confirm  this  hypothesis,  we  considered  the  next  case. 

2.  g^(x,y)  ^  (x  +  iy)/(l  +  X  +  iy).  This  case,  as  well 

as  the  theoretical  values,  was  suggested  by  0.  A.  Gross,  and 
he  determined  the  exact  values. 
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iiinction  |  Exact 

Computed 

1 

1.17 

1.17 

f2(1.0) 

0.667 

0.647 

8.  DISCUSSION 

As  can  be  seen,  the  agreement  in  general  is  quite  satis¬ 
factory.  We  can  obtain  reasonably  accurate  values  of  the 
maximum  return  and  of  the  optimal  allocation  policy  using 
small  amounts  of  machine  time. 

Combining  these  techniques  with  the  .method  of  the 
Lagrange  multiplier,  we  can  expect  to  solve  three- and  four- 
dimensional  resource  allocation  problems.  Extending  the 
method  to  cover  the  approximation  of  functions  of  5  and 

6  variables,  we  can  treat  Hitchcock— Koopmans  allocation 
processes  of  quite  high  dLmension. 

Finally,  if  we  combine  these  techniques — polynomial 
approximations  and  Lagrange  multipliers — with  that  of 
successive  approximations,  there  should  be  very  few 
allocation  processes  which  still  resist  our  efforts. 


u  • 
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Appendix  A 

1  OIHENSIONAL  ALLOCATION  VIA  DYNAMIC  PROGRAMMING. APPROXIMATIONS 

LIST  - 

•  LABEL 

C6A9SC1 

COMMON  N.M.JR.H.lH.L.Y.XL.Ii.A.P.GtNI.AI.Bl 
X  tlHl.Ll.LZ.Hl.S.Sl.HOLAST.HOLBL.HOLOOT 
CIMENSICN  XL(20).WI20),A(20),P(20).F|20).X{20I.A1I10).B1( 10) 

1  CALL  INPUT 
IF|N-100)2.3.3 

C  STAGE  1 

2  NI  >I 

DO  A  K*1.M 
SUM»0. 

DO  5  J«ltJR 
V^XLI J) 

CALL  RETG 
CALL  POLY 
FI J)*G 
XIJ)»Y 

S  SUM>SUM«hlJ)«F(J)*PIK) 

AlK)*SUM 
A  CONTINUE 

OUTPUT  A2t60tNI. JR.M 

60  FURMAT(20H1  ITERATIVE  STAGE  IS.5H  R-I3.5H  M«I3) 

OUTPUT  A2.66I 

661  F0RMAT(S3H0  DECISION  RETURN  AT  R  ROOT  VALUES  I 

OUTPUT  A2«6I.lX(J).FU)tJ«I»JR) 

61  FCRMATIIH  2E20.8) 

OUTPUT  A2.S50 

550  FORMAT  ( 30H0  COEFFICIENTS  A|K)  ) 

OUTPUT  A2.50.IA|K).K«1,M) 

50  FORMATIIH  E20.a) 

CALL  OUT 

C  stage  NI  GREATER  THAN  I 
DO  6  N2>2tN 
NI«N2 

DO  7  J*1,JR 

V=0, 

CALL  RETG 
Y»XHJ)-0. 

CALL  POLY 
SUP»0. 

CO  8  K-1,M 
8  SUP«SUM«AiK)«P|K) 

FIJ)»C*SUP 

X(J)«0. 

CO  9  I«2,IH 
XI»I 

IFI  |XI-I.)«H-A1IN2) ) 9, 110.110 

110  CONTINUE 

IF{(XI-I.)*H-GIIN2) )  111.  HI.  11 

111  CONTINUE 

IF((XI-1.)«H-XL(J))10.10.I1 
10  Y»|XI-1.)«H 
CALL  RETG 
Y«XLI J)-IXI-I.)*M 
CALL  POLY 
SUP»0. 

DO  12  K>1.M 
12  SUM>SUM«A(K)»P(K) 
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PHsG^SUI* 

IF(PH-FIJ) I9t9,l5 

15  FIJ)«PH 

9  CONTINUE 
11  CONTINUE 
7  CONTINUE 

OUTPUT  «2t60tNIt JRtH 
OUTPUT  A2,661 

OUTPUT  «2t61t(X( JI.FI J),J>1, JR) 

00  16  K*1(P 
SUP>0. 

00  17  J>1«JR 
VsXLiJ) 

CALL  POLY 

17  SUP«SUP«WfJ)«F(J)«P(K) 

AIK)«SUP 

16  CONTINUE 
OUTPUT  62tS50 

OUTPUT  62tS0t( AIKI .K^l.NI 
CALL  OUT 
6  CONTINUE 
GO  TO  1 
3  CALL  EXIT 
STOP 
ENC 

•  LIST 

•  LABEL 
C6A9S02 

C  INPUT  ALL  CONSTANTS 

SUBROUTINE  INPUT 

COPMON  NtPt JRtH,lH.LtY»XLtH»A,P,CtNIf AltBl 
X  t IHl>Ll.L2tHl(StSlfHOLAST»HOLBLtHOLOOT 
01  PENSION  XL(20).Nl20).At20).PI20),Fi20)tXl20).Alf lO).BI(lO) 
INPUT  Al,20«N.Pt JRt IH.L 

INPUT  41,21,|XLI J).J«1*JRI*(U| l)»I>l«JRI,H,(Alf lit |«1,N), 

X  (B14 J)t J-1»N) 

OUTPUT  A2«22«N.P,JR.1H,L 

OUTPUT  A2*23t(XL( J)tJ-ltJB)tlWlI),I«ltJRItCfH,|AI( ntI«ltN)» 
X  IB1IJ)»J«I.N) 

INPUT  61.1001, IHl, LI. L2 

INPUT  61. 1002. HI. S, SI 

INPUT  61.1000,HUt.AST.HOL6L,HOLC0T 

RETURN 

ICOO  F0RPATI3A1) 

ICOl  FORPATfllO) 

1C02  F0RPATIF10.3I 

20  FORPATISIIO) 

21  F0RPATIE20.8) 

22  FORMATIIHISIIO) 

23  FORPATIIH  E20.a) 

ENC 

•  LIST 

•  LABEL 
C669S03 

C  COPPUTE  RETURN  G  AS  FUNCTION  OF  V 

SUBROUTINE  RETG 

COPPON  N.P.JR.H.lH.L.V.XL.M.A.P.G.NI.Al.Bl 

X  , IH1,L1,L2.H|.S.S1.HOLAST.HOLBL.HOLOOT 
Cl  PENS  ION  XLl20).Nt20).AI20) .P(20).FI20).XI20).A1I lOI.BH 10) 


17- 


XNI^Nl 

C«XNI*SCRTF(Y>1.) 

RETURN 

ENC 

•  LIST 

•  LABEL 
C6A9S0A 

C  NORMALIEEC  POLYNOMIALS  PIKI  AS  FUNCTIONS  OF  Y 
SUBROUTINE  POLY 

COMMON  NtM. JRtHf IHtLtYtXLtMtA.PtCtNItAltBl 
X  ,IHltLl»L2,Hl,StSltHOLAST.HOLBLtHOLOOT 
DIMENSION  XLI20JtWl20).AI20I  •PI20)»FI20)#XI20I#AII  lOl.BHIOI 
P(1)>1. 

PI2l*SQRTFI3,)*(2.«Y-l.) 

Ml»M-2 

DO  40  I-ItMl 
BI>I 

0>SQRTFI2.*BI^3.)/(BI«1.) 

E«SQRTFI2.«8I«I.) 

B>BI/SQRTFI2.«8I>1.) 

40  Pn>2)»C»IE«U.»Y-l.)«PI  I»l»-B*PC  I)  » 

RETURN 

END 

•  LIST 

•  LABEL 
C649S0S 

C  COMPUTE  and  output  TOTAL  RETURN  FUNCTION  AT  STAGE  NI 
SUBROUTINE  OUT 

COMMON  NtMt JRtHt IHtLtYtXLtWtA*P«GtNItAlt61 
X  «IHl,Ll,L2tHI,StSltHOLASTtHOLBLtHOLOOT 
DIMENSION  XLI20) tNl20) f AI20) tPI20>f FI20)#Xl20If AH lOIfBlI lOI 
DIMENSION  ZIIIO) 

OUTPUT  42*7StNlt JRtM 

75  F0RMATI20HI  ITERATIVE  STAGE  15, 5H  R-I3,5H  N-I3) 

IF(L-1)13,18,18 
13  OUTPUT  42,774 

774  FORMAT  1 40H0  RESOURCE  RETURN 

DO  72  IS»l,IH 
XI»I9 

Y»«XI-1. )»M 
CALL  POLY 
SUM«0. 

DO  73  K»l,M 

73  SUM«SUM«AtK)*P|K) 

FF*SUM 

OUTPUT  42,74,Y,FF 

74  FORMATIIH  FI5.3,E20.8I 
72  CONTINUE 

GO  TO  71 
18  CONTINUE 

CO  772  I9>1,IH1 
XI-I9 

Y»!XI-l.l*Ml 
CALL  POLY 
SUM«0. 

DO  773  X-1,M 
773  SUM»SUM^A<K»*PIX) 

FF»SUM 

N5*FF«S*Sl 

no  1003  I5>1.110 
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ICC3  Z(I5)=bCLBL 

CG  1004  L5=1,L2 

[F( !9-Ll*L5-ll?6,77, 76 

77  DO  78  16  =  1 ,  1 10 
Z ( 16 ) =HCLnCT 

78  CCNTINUE 
76  CGMINUb 

1C04  CONTINUE 

DC  1005  15=1,110,10 
1CC5  Z ( 15) =HCLCCT 
Z (N5)=HCLAST 

OUTPUT  42, 1CC6, Y,  ( Z (  1 7) ,  i 7= 1 , 105  ) 
772  CONTINUE 
71  CONTINUE 
RETURN 

1C06  FORMAT!  IH  F  10.3, ICOAl) 

FNC 
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Appendlx  B 

C  2  DIMENSIONAL  ALLOCATION  VIA  DYNAMIC  PROGRAMMI NG» APPROX IMAT IONS 

*  LIST 

*  LABEL 
C649500 

COMMON  NtM.JR.NI tX . Y .G«F 1 .W* V»XL t YL *A *P *PH«P«B*C • IH tH. A1 *B1 (C 1 *01 
DIMENSION  W(20) tV(20) *XLI20) tYL(20) •A(20*20)«P(20)«PH(20I •FI20*20I 
X  tBI3) >C(3) •SUM(20> •SUM2(20) tPX(20t20) •PY(20*20I 
X  •ai(20)tAl(20)»Cl(20)*01(20) 

1  CALL  INPUT 
IF(N-100J2*3.3 

C  STAGE  1 

2  Nisi 

OUTPUT  42*60*NI (JR»M 

60  F0RMAT(20H1  iTERATiVE  STAGE  I5.5H  RsI3»5H  MsI3) 

OUTPUT  42*54 

54  FORMAT (90H0  X  Y  F 

X  XDEC  YDEC  I 

DO  4  1=1* JR 
XsXL(  I ) 

call  poly 

DO  600  K«ltM 

600  PX(K.I)=P<K) 

DO  41  JsltJR 
YxYL( J) 

CALL  PHOLY 
DO  601  L»1.M 

601  PY(L.J)=PHIL) 
call  RETG 

F( 1 ,J)=G 
XDECsXLI I ) 

YDECsYL(J) 

OUTPUT  42.61,X.Y.Fn.J»  .XDEC.YOEC 

61  FORMATdH  3E20»6 .2F 10. 3  ) 

41  CONTINUE 

4  CONTINUE 
DO  5  K=1»M 
DO  51  L=1.M 
DO  52  1=1. JR 
SUM( I )=C, 

DO  402  J=1,JR 

402  SUM(  I  )=SUM(  1  )-*V(  J)*F(  I  ,J)*PV»L.J) 

52  CONTINUE 
SUM1=0. 

DO  7:  1=1, JR 

70  SUMl=SUMl*-w(  I  )*PX(<,I  )*SUM<  I  ) 

A(K,L)=SUM1 
51  CONTINUE 

5  CONTINUE 
OUTPUT  42.50 

50  FORMATOOHl  COEFFICIENTS  A(K,L)  ) 

OUTPUT  42.53. ( « A ( I ,J) ,I»1,M) . J»1 .M) 

53  FORMATdH  20.6) 

OUTPUT  42.556 
DO  55  1=1.3 

X  =  bd) 


-eo- 


V*C«I) 

CALL  OUT 
55  CONTINUE 

C  STAGE  N1  GREATER  THAN  1 
DO  6  N2>2*N 
NI>NI^1 

OUTPUT  42*60*Nl«JRiM 
OUTPUT  42*54 
00  6  Isl,jR 
DO  80  J«1*JR 

x»o.  \ 

Y»0.  ^ 

call  retg  \ 

X«XL«I)-0. 

Y«YH  JI-0. 

CALL  TOTRET 
Ft  I f JJaFl 
X0EC>0* 

YDEC*0. 

00  9  IlsUIH 
X1>I1 

IFKXI-I.  )*H-AU  n  >9*99*99 
99  CONTINUE 

IFr  «X1-1.)*H-B1( I) >991*991*11 

991  CONTINUE 

IF( |XI-1.>»H-XLI I ) >10*10*11 

10  CONTINUE 

DO  91  JlaUlH 
XJ»J1 

IFHXJ-1.>»H-CU  J>  >91*992*992 

992  CONTINUE 

IF< (XJ-1*>*H-D1(J) >993*993*9 

993  CONTINUE 

IFJ IXJ-1. >*H-YH J) >100*100*9 
100  X=IXI-1.>*H 
Ya(XJ-l.>»H 

call  retg 

X»XLJ I >-(Xl-l*>»H 
Y*YL( J>-(XJ-1. >*H 
CALL  TOTRET 
F2»F1 

IF|F2-F( I ,J>  >9: *91*15 
15  F( I,J>*F2 

X0ECx(Xl-l.>*H 
YDECx(XJ-l.)*H 
91  CONTINUE 
9  CONTINUE 

11  CONTINUE 

OUTPUT  42*61*  XL( I > *YL( J) *Fn .J> *XOEC*YDEC 
80  CONTINUE 
8  CONTINUE 
DO  400  Xal*M 
DO  401  L«1*N 
DO  502  I>1*JR 
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SUM« I )=0, 

00  500  Jsl.JR 

500  SUMJ I )=SUM( I )fV( J»*F( I ♦J)*PY(L.J) 

502  CONTINUE 

SUMlsO. 

DO  501  1=1 .JR 

501  SUMl  =  SUMl-*-WI  I  )*SUM<  1  )*PX(Ktn 
A(IC*L)=SUM1 

401  CONTINUE 
400  CONTINUE 

OUTPUT  42.50 

OUTPUT  42.53. ( (A( I »J) »I*1.M) ,J»1»M) 

OUTPUT  42.556 

556  FORMAT (54H0  TEST  X  Y  F  > 

DO  555  1=1.3 
X  =  B( I  I 
Y  =  C(  I  ) 

CALL  OUT 
555  CONTINUE 
6  CONTINUE 
GO  TO  1 
3  call  exit 
STOP 
END 

*  LIST 

»  LABEL 

C649501 

C  INPUT  ALL  CONSTANTS 

SUBROUTINE  INPUT 

COMMON  N.MtJR.M .X.Y.G.FI tW.V.XL.YL.A.P.PH.F.B.C.IH.H.Al .Bl.Cl.Dl 
DIMENSION  W(20) .V(20) .XL (20) .YL(20) •A(20.20) .P(20) .PH(20  > .F(20.20) 
X  .B(3) .C(3) .SUM(20) .SUM2(20).PX(20.20) .PY( 20.20) 

X  .B1(2C).A1(20).C1(20). 01(20) 

INPUT  41.20.N.M.JR.IH 

INPUT  41 .21.(XL( J) .J=l .JR) .( W( I ) .1  =  1 .JR ) . ( YL( J) .J=l  .JR)  . 

X  (V( 1 ) . I>1.JR) .H.IBI I ) .I=1.3).(C( J) .J«1.3) 

OUTPUT42.22.N.F  .JR.IH 

OUTPUT42.23. (XL( J) . J- 1 . JR ) . ( W ( I ) . I = 1 . JR ) . ( B ( I ) . 1 « 1 . 3 ) .H 
INPUT  41. 233. (Aid). 1  =  1.  JR). (Bid). 1=1.  JR  ).(C1(I  )  .1  -  l.JR)  . 

X  (Old). 1  =  1. JR) 

233  FORMAT(F10.3) 

20  F0RMAT(4I1C) 

21  FORMAT(E20.8) 

22  FORMAT) 1H141 10) 

23  FORMATdH  E20.8) 

RETURN 

END 

*  LIST 

*  LABEL 
C649502 

C  COMPUTE  RETURN  G  AS  FUNCTION  OF  X  AND  Y 

subroutine  RETG 

COMMON  N.M.JR.MtX.Y.G.Fl.W.V.XL.YL.A.P.PMtF.B.C.IM.M.Al.Bl.Cl.Ol 
DIMENSION  w(20;.V(20).XL(20).YLt20).A(20.20).P(20).PM(20).F(20.20) 
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X  *B(3)  (CO)  tSUH(20)  .SUM2  ( 20  )  tPX  (  20 120  )  tPY  (  2C  ^20  I 

X  tBl(20)tAl(20)tCl(20)t01(20) 

XNI=NI 

GaSQRTF(2.*XNI-l.)*SQRTF(SQRTFl ( X+. 1 )* ( Y^. 1 ) »  ) 

RETURN 

END 

*  LIST 

*  LABEL 
C649503 

C  NORMALIZED  POLYNOMIALS  P(<)AS  FUNCTIONS  OF  X 

subroutine  poly 

COMMON  N.MtJR*NI >X » Y.G.F 1 *W*V*XL t YL .A.P *PHtF .BtC » IH tH t A1 »B1 tCltOl 
DIMENSION  W(20l *V(20) •XL<20>  >YL(20) tAIZOtZO) iPIZO) tPH(20 I (FI  20*20) 
X  »B(3)  .0  3)  »SUM(20)  .SUM2(20)  .PX(20.20>  .PYI20.20) 

X  .Bl (20) •A1(2C) *C1 (20) .01 (20) 

P( 1  )  =  1. 

P(2)aSQRTF(3, )*(2.*X-1.) 

MlaM-2 

DO  40  1 1=1. Ml 
Blall 

OaSQRTF(2.*BI-*-3.  )/(BI  +  l.  ) 

E»SQRTF(2.*dI-.l.  ) 

TaBI/SQRTF(2.*£,I-l.) 

40  P(  1 1^2)=D*(E*(2.*X-1.  )»P(  1 1-f  1  )-T*P(  II )  ) 

RETURN 

END 

*  LIST 

*  LABEL 
C649S04 

c  normalized  polynomials  ph(X)  as  Functions  of  y 

SUBROUTINE  PHOl Y 

COMMON  N.M.JR.NI *X . Y *G*F 1 , W* V  *  XL  * YL * A *P *PH*F .B.C *  I H .H* A1 .Bl.Cl.Dl 
DIMENSION  W(20) .V(20) .XL(20).YL(20) .A(2C.2C) .P(2C) .PH(20) .F(20.20) 
X  .B(3)  .0  3)  .SUM(20)  .SUM2(20)  .PX(20.2C)  .PY(20.20) 

X  *B1(20).A1(20)*C1(23)*01(20) 

PH(1)=1, 

PH(2)=SaRTF(3,)*(2.*Y-l. ) 

Ml=M-2 

DO  90  Ilal.Ml 
31*11 

D=SORTF(2.*BI+3. )/(BI>l. ) 

E  =  SQRTF(2.*ai'H.) 

T=BI/SQRTF(2.*BI-1.) 

90  PH(  I  l  +  2)aD*(E*(2.»Y-l.)*PH(  11-.1)-T«PH(  ID) 

RETURN 

END 

*  LIST 

*  LABEL 
C649S05 

C  COMPUTE  total  return  AS  FUNCTION  OF  X.Y 

subroutine  TOTRET 

COMMON  N.M.JR.NI *X . Y *G.F 1 . W . V .XL  * YL . A *P .Ph *F *B *C *  1 H .H . A1 .Bl.Cl.Dl 
DIMENSION  W(20) *V(20>  *XL(20) .YL(20) .A ( 20 *20 ) *P ( 20 ) *PH ( 20 ) .F) 20*20) 
X  .B( 3) .O  3 ) .SUM (20) .SUM2 (20) .PX( 20.2  0) .PY( 20.20) 
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X  .Bll20)*Al(20).Cl<20)tDl(20) 
CALL  POLY 
call  PHOLY 
DO  81  IL*1*M 
SUM2( IL|30. 

DO  83 

83  SUM2<  IL)=SUM2(  IL  J^P  (  lit » *A (  IK »  IL  ) 
81  CONTINUE 

SUM3=0. 

DO  8A  ILsl*M 

84  SUM3  =  SUM3-*’PH(  IL)*SUM2(  IL) 
F1sG4'SUM3 

RETURN 

END 

*  LIST 

*  LABEL 


C649506 


TOTAL  RETURN  FROM  POLYNOMIAL  APPROX  I  MAT  I  ON tFUNCT I  ON  OF 
SUBROUTINE  CUT 

COMMON  N»M»JR»M»X#Y»G«FltW#V»XL«YL#AtP»PHtF»B»CtIH»Hi 

DIMENSION  W«20)»V(20) .XL ( 20 ) .YL t 20 ) .A t 20 .20 ) »P ( 20 ) »PH( 
X  ,B(3) .C(3I .SUMI20) .SUM2t20) .PX(20.20) .PY(20»20) 

X  ,B1(20).A1(20).C1(20).DH20) 


X.Y 


Al.Bl.Cl.Dl 

20)»FI20»20) 


CALL  POLY 
CALL  PHOLY 
DO  281  IL*1.M 
SUM2( IL>*0. 

DO  283  IK«1.M 

283  SUM2(  IL)*SUM2(  IL  l■♦•P(  IK)«A(  IK.IL) 
281  CONTINUE 

SUM3=0. 

DO  284  IL«1.M 

284  SUM3»SUM3>PH( IL)*SUM2( ID 
FF*SUM3 

OUTPUT  42.261 .X.Y. FF 


RETURN 

261  FORMAT! IH  3E20.6) 
END 
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